Virial theorems for vortex states in a confined Bose-Einstein condensate 
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We derive a class of virial theorems which provide stringent tests of both analytical and numerical 
calculations of vortex states in a confined Bose-Einstein condensate. In the special case of harmonic 
confinement we arrive at the somewhat surprising conclusion that the linear moments of the particle 
density, as well as the linear momentum, must vanish even in the presence of off-center vortices which 
lack axial or reflection symmetry. Illustrations are provided by some analytical results in the limit of 
a dilute gas, and by a numerical calculation of a class of single and double vortices at intermediate 
couplings. The effect of anharmonic confinement is also discussed. 

PACS numbers: 03.75.Lm, 47.32. Cc, 47.37.+q 



Quantized vortices observed in a bulk superfluid such 
as 4 He have fascinated physicists for a long time |]| be- 
cause they provide definite macroscopic manifestations of 
subtle quantum phenomena. The subject has been sig- 
nificantly enriched in recent years with the realization of 
ultracold atomic Bose-Einstein condensates (BECs) con- 
fined in a finite region, where the strength of effective 
interactions may be manipulated by varying the number 
of atoms in a given trap. Typically, these condensates 
are sufficiently dilute that a mean-field Gross-Pitaevskii 
(GP) approximation is reliable. It is then possible to 
explicitly calculate a variety of vortex states which are 
relevant to experiment. There have been numerous con- 
tributions in this area, partly reviewed in Ref. [2j, but 
the subject is still active because a number of finer issues 
remain unexplored. 

The GP approximation is adopted throughout this pa- 
per. We mainly consider an effectively two-dimensional 
(2D) Bose gas of N atoms, each with mass M, which 
interact pairwise with a contact potential of positive 
strength Uq and are confined by an axially symmetric 
external potential tkuoV(p/ao). Here the constant loq 
carries dimensions of frequency, a — \/?i/Mujq is the 
corresponding oscillator length, and p — \J x 1 + y 2 is the 
radial distance from the center of the trap. Rationalized 
units are introduced by measuring time t in units of 1/ujq, 
distances x and y in units of clq, while the condensate 
wave function is rescaled according to 'F — * VN <F/ao 
and thus acquires unit norm: J dxdy = 1. The en- 
ergy functional is then given by 



E= I - (V>F*V<F) + ^( i o)<F*'F + ^(<F*>F) 2 dxdy, 

and yields energy in units of HloqN; V{p) is the rational- 
ized trap potential and g = MNUo/h is a dimensionless 
coupling constant. 

In a frame rotating about the center of the trap with 
constant angular frequency u> (in units of ujq) stationary 



states of the gas satisfy the time-independent differential 
equations 



/i^ — lUJEapXadp^ — 

/i'F* + iuj£ a px a d(3^* = 



SE 
Sty* 
SE 



where 



SE 



= — A* 
2 



V(p)m + g 



(2) 
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and 5E/5^> is its complex conjugate. A chemical poten- 
tial jU (in units of huo) is introduced in Eq. J5J in order to 
enforce a definite number of particles. Greek indices a, (3 
take two distinct values corresponding to the two spatial 
coordinates x\ = x and x-i = y, and the Einstein summa- 
tion convention of the repeated (dummy) indices is con- 
sistently employed throughout the paper. Finally, e a /3 is 
the usual 2D antisymmetric tensor and d a — d/dx a . 

Now, given a solution 'F = "F(x, y\p, uS) of Eqs. J2J), the 
time dependent wave function 

*{x,y,t) = ^(x^y'l^^e-^, (4) 
x — x cos Lot + ysinwi, y — — x sinwt + y cos tot, 

satisfies the standard Gross-Pitaevskii equation in the 
laboratory frame (id^/dt — SE/S^>*) and may be 
thought of as a configuration that rotates (precesses) 
about the center of the trap with angular frequency ui. 
While all calculations will be based on the stationary 
Eqs. J2J, Eq. Q is important for a proper interpretation 
of the results. 

Before discussing explicit solutions, we derive a class 
of virial theorems which follow directly from Eqs. J2J- 
Thus we multiply both sides of the first equation by , 
the second by 'F, and add the two equations. We then 
integrate both sides over the entire xy plane and apply 
partial integration to obtain 



p + uj£ = E'ki, 



E; 



trap 



2E, 



pot: 



(5) 
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where -Ekin, -Strap and E pot correspond to the three terms 
in the total energy E of Eq. In the left hand side of 
Eq. 10) we have employed the relations 

J dxdy = 1, I = - J ^*e a px a dp^ dxdy, (6) 

where the first is consistent with our choice of rational- 
ized units and the second is the definition of the angular 
momentum per particle (in units of K). Eq. (J5J is the 
first of a series of virial relations that must be satisfied 
by all solutions of Eqs. (J2J. 

We now repeat the procedure by multiplying the first 
Eq. (J2J by the second by and then subtracting the 
two equations to obtain 



uJE al 3X a dpn + d a J a = 0, 



J 



(7) 
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where n and J are the familiar particle and current 
densities. The same result could be derived by ap- 
plying the continuity equation in the laboratory frame 
(dn/dt + V • J = 0) for a wave function of the form 10}. 
Now multiply both sides of Eq. JJJ) by x a and then apply 
partial integration to obtain the virial relation 



P a + LOSapRp = 0, 



(8) 



J a dxdy, R a 



x a n dxdy, 



where P = {P\,P2) is the linear momentum in the ro- 
tating frame and R = (Ri, R2) may be thought of as the 
mean position of the configuration in question. It should 
be clear that neither P nor R is conserved in the labo- 
ratory frame, where they precess about the center with 
frequency w, in complete analogy with the momentum 
and position of a pointlike particle in circular motion. In 
view of this analogy the virial relation JSJ) appears to be 
quite natural. 

A third and more elaborate class of virial relations is 
obtained starting again from Eqs. J2J) but now multiply 
the first equation by d a ^!*, the second by d a ^f, and then 
add the two equations to yield after some rearrangement 



8E 5E 
pd a n + ux al = — <9 Q ** + — d a ^ = f a (9) 

where n = '3>*\I> is again the particle density, 



1 



7 = -e a0 d a ^*d^ 



(10) 



may be referred to as the topological vorticity 0- [3- 13 ■ 
and we further use Eq. (|3J) to write 

fa = dp<7 af 3 + V(p)d a n, 

o aP = w5 a p-^(d a y*d^ + d **d a V) : (11) 

w = \ [(V** ■ V*) +. 9 (1'**) 2 ] . 



Hence Eq. © reduces to a more transparent form: 
fid a n + u x a j — djia a fi -I- V(p)d a n, 



(12) 



which will provide the basis for the derivation of a num- 
ber of interesting virial relations. 

Some key elements of the preceding discussion, such as 
the topological vorticity 7 and the tensor <j a a, appeared 
earlier in a study of the magnetic continuum ylOl as well 
as of homogeneous superfluids 0]. In the latter case, the 
total topological vorticity r = J 7 dxdy = 2im is integer 
valued (n — 0, ±1, ±2, . . .) for wave functions that satisfy 
the boundary condition |^| — > 1 at spatial infinity. How- 
ever, for the confined gas under present consideration, 
the relevant wave functions satisfy the boundary condi- 
tion I 'J I — * which leads to Y = by a straightforward 
partial integration. Similarly, the linear and angular mo- 
menta defined from 



e af3 xp^dxdy, I 
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p 2r ydxdy, (13) 



may be shown to coincide with the standard definitions 
given in Eqs. JSJl and JSJl by freely performing partial in- 
tegrations which are fully justified in a confined gas. In 
contrast, partial integrations are generally ambiguous in 
a homogeneous gas and Eqs. (|T3)l do not coincide with 
the standard definitions, thus leading to a subtle distinc- 
tion between momentum and impulse 0, Q . Certainly, 
for our current purposes, Eqs. (|13|l may be employed in 
conjunction with Eq. (|12|l without further questioning. 

Thus we integrate both sides of Eq. (|12fl over the en- 
tire xy plane and note that the terms d a n and dpaap 
lead to vanishing surface integrals at spatial infinity. The 
remaining terms may be arranged to yield the virial re- 
lation 



P a + £ a —V'n dxdy = 0, 
J P 



(14) 



where we have employed the linear momentum P a from 
Eq. (|13f) and performed a partial integration in the second 
term, with V = dV/dp. 

The basic relation (|12f> may be further iterated by mul- 
tiplying both sides with xp and then integrating over all 
space to obtain 



p8 a [3 — Lo J x a xp^dxdy = 

J <r a p dxdy + J (^V8 a p H — ^-J^-V'^j n dxdy, (15) 

which may be applied for any combination of indices a 
and (3 and thus contains three independent virial rela- 
tions. An interesting special case is obtained by taking 
the trace of both sides of Eq. (flTf : 



2p — lo J p 2 ^ dxdy = 

J tradxdy + J (2V + pV) ndxdy, (16) 
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where we may further insert the definition of the angular 
momentum I from Eq. 1|13[1 and tra = g(^*^) 2 from 
Eq. ifrrjl to write 

p + tu£ = Spot + f (V + ~ p V) n dxdy. (17) 

This virial relation may be derived also by applying a 
Derrick-like Q scaling argument to the extended energy 
functional F = E — ujL — pN. 

Eqs. ©, (jSJl, {II} and {T7|l already provide an inter- 
esting variety of virial theorems which are employed in 
the following to check and analyze explicit solutions of 
Eqs. It should be noted that virial relations 101 and 
© are insensitive to the specific choice of the trap poten- 
tial, while l|14f) and l|17|l depend crucially on the choice 
of V = V(p). Similarly, Eqs. JHJ and (|14|l are insensitive 
to the specific form of the tensor a a p and are thus valid 
for any type of interparticle interactions. 

Most of the theoretical models employed to describe 
realistic BECs assume a harmonic trap potential V — 
\p 2 and, hence, V = p. The virial relation (fT7|) may 
then be written in the form p + uj£ — S pot + 2S tra p which 
is combined with Eq. (JSJ) to yield 

Skin + Sp ot = S trap . (18) 

This relation does not contradict the existence of finite- 
energy stationary solutions in the rotating frame and 
must indeed be verified by any such solution of Eqs. 
Similarly, the virial relation 114JI simplifies for V' = p 
to read cuP a + e a pRp — which may be combined with 
Eq. JSJ) to arrive at the somewhat surprising conclusion 
that both the linear momentum and the linear moments 
of the particle density must vanish in a harmonic trap: 

P a = = R a , (19) 

provided that u> ^ 1, a restriction that is not essential 
in the case of repulsive interactions because stationary 
states are then possible in a harmonic trap only for u < 
1(= u>o). In the special limit lu — 1, which could be 
achieved in the case of attractive interactions Q, P and 
R need not vanish. 

A simple explanation of the preceding result can be 
obtained by noting that for a system of atoms in a har- 
monic trap with translationally-invariant interparticle in- 
teractions, the center-of-mass (CM) coordinate separates 
from the internal coordinates and behaves as a free par- 
ticle in a harmonic well. For to < 1, the CM must be 
in its ground state, so P and R must be zero because 
they depend only on the CM coordinates. At w = 1, 
it is possible to put the CM into a rotating state with 
nonzero P and R. This is relevant for the case of attrac- 
tive interactions where the rotating states with lu = 1 
are such that the angular momentum is carried by the 
center-of-mass This viewpoint makes it clear that 
the same result will apply for any translationally invari- 
ant interparticle interactions, in the case of 3D harmonic 
confinement, and beyond the CP approximation (for the 



expectation value of the CM and conjugate momentum). 
It is also clear that Eq. (19) is not valid in the case of 
anharmonic confinement, as discussed later in this paper. 

Explicit solutions of the Gross-Pitaevskii theory were 
initially obtained in the limit of a very dilute gas |lOj| . 
A wave function with definite angular momentum £ may 
then be restricted to the lowest Landau level (LLL): 

E z rn e -\z\ 2 /2 
c m * m , * m = , (20) 

m>0 

where the sum extends over nonnegative integer m (for 
positive £) and z = x + iy. The unknown coefficients 
c,„ are calculated by minimizing the total energy S = 
£ + g/2 J ('I'**) 2 dxdy under the constraints yi \c m \ 2 = 1 
and J2m\c m \ 2 = I. This task was initially |l0j| carried 
out numerically to furnish an impressive variety of vor- 
tex states leading up to a vortex lattice for lar ge £ . For 
small £, some analytical results were obtained [ll| by a 
perturbative expansion of the coefficients c m in powers 
of £ or £ = 1 - I. 

These results may already be used to illustrate the 
virial relations l|19|l . We first consider the linear moments 
of the particle density R a , with a — 1 or 2, or their com- 
plex combination Ri+iR 2 = J *&*z i S> dxdy. We may then 
insert the series representation (|20|l for the wave function 
^, note that z^ m — \Jm + l\& m+ i, and apply the usual 
orthogonality relations for the ^m's to obtain 

i?i + iR 2 = ^2 VrnTTc* m+l c m . (21) 

m>0 

If we now use the perturbative expansions for the c m 's 
from Ref . ^lj ; with due attention to phase (sign) conven- 
tions 01, we find that R\ + iR 2 — order-by-order in a 
consistent expansion in powers of £ or £ = l—l. A similar 
calculation of the linear momentum yields P\ + iP 2 = 0, 
thus confirming the validity of both virial relations in 
Eq. (|19fl , as well as providing a nontrivial check of con- 
sistency of the results of Ref. 

El 

A more convincing demonstration is possible over the 
entire range < t < 1 where a closed-form expression 
for the optimal LLL wave function was recently achieved 

m 

pi 

y = il u x - b ) + iy] e ~h[(*~*f-^y+y 2 ] 

a=(V~£-£) 1/2 , b=— -, (22) 

a 

which describes an off-center vortex located on the x axis, 
modulo an overall azimuthal rotation, at a distance b = 
b(£) from the center of the trap. For £ — > (b — > oo) the 
vortex is expelled from the system, while for £ — > 1 (b — > 
0) the vortex moves to the center and becomes axially 
symmetric. The corresponding particle density reads 

n = ^ll( x - b )2 + y2] e -(x-") 2 ~y 2 (23) 
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and lacks axial or reflection symmetry except for £ — 
or 1. Nevertheless, an explicit calculation shows that 

Rl = J xndxdy = V£[2a-b + a(a-b) 2 }=0, (24) 

using the explicit expressions for a and b from Eq. (|22[) . 
and i?2 = J yndxdy = 0, thanks to the y — > — y symme- 
try of the particle density. A similar calculation shows 
that the linear momentum P = (-Pi,/^) also vanishes, 
thus verifying both virial relations in Eq. (|19|) in spite 
of the lack of axial or reflection symmetry in the wave 
function (1221) . This result appears to be surprising if we 
naively view the off-center vortex as a pointlike particle 
rotating on a circle with radius b. In fact, the spatial dis- 
tribution of particle and current densities is more subtle 
in an off-center vortex and leads to vanishing P and R, 
even though the vortex does precess about the center of 
the trap. As we shall see shortly, this curious result is 
valid only for harmonic confinement (V = \p 1 \ 

We complete the discussion of the dilute-gas limit by 
quoting the energy E and frequency u) associated with 
the wave function l|22() : 

with S = g/87T, a result that is valid in the limit <5 <C 
1 . A notable feature of this limit is that frequency uj is 
independent of angular momentum I [jj, ITflL ITU H^ j . 

One should keep in mind that practically all experi- 
ments have been performed on BECs with S > 1, of- 
ten S 3> 1, where the weak coupling (LLL) theory is no 
longer valid. A numerical solution of Eqs. (J2J) is neces- 
sar y fo r strong couplings. A numerical method developed 
in [14| is based on a norm-preserving relaxation algo- 
rithm which, in effect, capitalizes on the virial relation 
© to find a wave function of unit norm that is a (local) 
minimum of the energy functional in the rotating frame: 
E ro t = E — ui£. Here we employ a variation of the norm- 
preserving algorithm to minimize instead the Lyapunov 
functional E' = E + \a{£ — b) 2 where a and b are arbi- 
trary constants with a > 0. Local minima of E' satisfy 
Eqs. J5J with frequency determined self-consistently from 
uj = a(b — £) and chemical potential \x from Eq. ©. The 
constants a and b are chosen (tuned) to ensure conver- 
gence to nontrivial solutions with angular momentum £ 
in the desired range. The advantage of this algorithm is 
that it finds solutions of Eqs. J2J which are stationary 
points but not necessarily local minima of the functional 
E ro t ■ 

In the following we describe a class of solutions for 
the intermediate coupling S = g/Sn = 2 where the LLL 
approximation is quantitatively inaccurate. Needless to 
say, all virial relations, including Eqs. (|18|l and (|19|) were 
confirmed by our solutions to within numerical accuracy. 
In Fig. n we present contour plots of the particle density 
n for four characteristic values of angular momentum in 
the range < £ < 2, whereas Fig. [21 shows the results for 
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FIG. 1: Contour plots of the particle density n for a class 
of vortex states calculated in a harmonic trap with 8 = 2. 
The four characteristic values of the angular momentum I 
correspond to the four branches AB, B'T, TA, and A... in 
Fig. H 




FIG. 2: Angular frequency u (in units of u>o) as a function of 
angular momentum per particle I (in units of K) for a class of 
vortex states calculated in a harmonic trap with 5 = 2. 



the frequency dispersion ui = dE/d£ = oj{£) throughout 
the same range. 

For < £ < 1 the calculated configuration is an off- 
center vortex with energy E(£) that is a concave func- 
tion of £. Thus the frequency u> = dE/d£ is a decreasing 
function of angular momentum taking values in the fi- 
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nite range lub < w < ^>A with lob = 0.44 for i — > 1~ and 
= 0.62 for I — ► 0. The vortex precesses faster the far- 
ther it is located form the center of the trap {£ — > 0). In 
the opposite limit (■£ — > 1 _ ) the vortex moves to the cen- 
ter and becomes axially symmetric. There is no sense of 
precession in such a vortex because the limiting frequency 
lub may then be absorbed into an effective chemical po- 
tential p, = fi + ujb and the wave function Q reduces 
to a quasi-static configuration with chemical potential p,. 
Also note that the band of allowed frequencies reduces 
to a single point uj a = Wb = 1 — 5 in the dilute-gas limit 
(5 <C 1), as is evident from Eq. 

As the angular momentum increases beyond unity the 
energy remains continuous, but its first derivative ex- 
hibits a finite jump which leads to a new limiting fre- 
quency ujb' = 0.74 for £ — > 1 + . The original vortex 
becomes again displaced from the center for I = 1 + and 
a second off-center vortex appears at an asymmetric po- 
sition on the opposite side of the trap. This picture re- 
mains largely correct in the region 1 < I < 1.55 and 
leads to branch B'T in Fig. [21 with limiting frequencies 
ujb' = 0.74 and ujy — 0.66. Energy is again concave in 
this region and thus frequency is a decreasing function of 
angular momentum. 

At point r {£ = 1.55) the calculated configuration be- 
comes a reflection-symmetric two- vortex state where the 
two vortices are located at the same distance on opposite 
sides from the center of the trap. Such a symmetric state 
persists throughout the branch TA (1.55 < £ < 1.78) 
with corresponding frequencies in the range 0.66 < lu < 
0.76. Incidentally TA is the only branch where the en- 
ergy E{£) is convex and thus the frequency u)(£) is an 
increasing function of angular momentum. 

Beyond point A (I > 1.78) the frequency becomes once 
again a decreasing function of angular momentum. This 
region seems to be characterized by the appearance of a 
new pair of vortices symmetrically displayed along the y 
axis, as indicated by the fourth (£=1.98) entry of Fig.^ 
Nevertheless, reflection symmetry appears to persist in 
this region. A related interesting question is whether 
or not hysteresis sets in when we reverse the cycle by 
reducing the angular momentum from, say, t — 1.98. 
While the cycle is perfectly reproduced down to point T, 
a reflection symmetric two- vortex state persists for some 
range of angular momenta below £r = 1.55, as indicated 
by the dashed line in Fig. [21 

The issue of stability of the calculated vortex states is 
rather delicate and may well depend on the specific ex- 
perimental protocol. According to ^} mechanical sta- 
bility requires that the energy be a convex function of 
angular momentum: duj/d£ = d 2 E/d£ 2 > 0. This condi- 
tion is satisfied only by the TA branch of Fig. [21 which 
corresponds to symmetric two-vortex states. In particu- 
lar, the whole of the AB branch, which corresponds to 
single off-center vortices, does not satisfy the criterion of 
mechanical stability. Nevertheless, precessing off-center 
vortices have been observed experimentally in a spheri- 
cal trap ^5|- Although the current two-dimensional cal- 
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FIG. 3: Angular frequency uj as a function of angular momen- 
tum per particle £ for a class of vortex states calculated in an 
anharmonic trap with 5 = 2 and A = 1 /4. 

culation does not directly apply to a spherical trap, a 
similar three-dimensional calculation [l6( leads to a class 
of U-shaped off-center vortices whose frequency disper- 
sion is completely analogous to the AB branch of Fig. [21 
Furthermore, the frequencies of precession measured in 
the experiment lie within the calculated frequency band 
[wb,o->a]- To conclude this digression, we note that the 
virial relations (jl9(l are valid also in a three-dimensional 
axially symmetric harmonic trap, with a = 1 or 2 corre- 
sponding to the directions perpendicular to the symme- 
try (rotation) axis, to be completed with P3 = = R3 
along the same axis. 

Finally, we consider the effect of anharmonic confine- 
ment modelled here by the rationalized trap potential 

V= ip 2 (l + A P 2 ), V' = P + 2\p\ (26) 

which is thought to describe the trap used in the ex- 
periment of Ref. 01 with A ~ 10~ 3 . In our numerical 
calculation we adopted a much larger A in order to em- 
phasize some generic features of anharmonicity. We have 
thus repeated our earlier calculation of vortex states in 
a harmonic trap (A = 0) now for A = 1/4 but the same 
coupling constant 5 = g/Sir = 2. 

The calculated frequency dispersion is shown in Fig.[3| 
which differs from Fig. [21 mainly by the fact that the dis- 
persion now extends well beyond u) = 1(= luq) because 
the condition lu < 1 is no longer necessary to ensure 
stability of the rotating gas. At first sight, the vortex 
configurations that correspond to the various branches 
of Fig. [3 are also similar to those calculated for the har- 
monic trap and shown in Fig. ^ However, a closer look 
reveals some subtle differences which are best illustrated 
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FIG. 4: Linear momentum P (in units of h/ao), linear mo- 
ment R (in units of ao), and generalized moment Q (in units 
of ag), as functions of angular momentum per particle I (in 
units of h), for single off-center vortices in an anharmonic trap 
with 5 — 2 and A = 1/4. 



by recalling the virial relation (|14(l now applied for the 
potential V of Eq. gg) : 



wP a + e a pRp = —2\e a pQf3, 
x a p 2 n dxdy, 



(27) 



where Q a is a generalized (higher) moment of the particle 
density. We further recall virial relation which is 
valid for any trap potential, and combine it with Eq. |2Tj) 
to yield 

Pa = n ZafiQfl, Ra = ~~ o Qa, (28) 

1 — W* 1 — W Z 

which differ significantly from Eqs. (|19|) in that the linear 
momentum P and moment R are no longer forced to 
vanish. 

For a numerical illustration we consider the class of 
single off-center vortices, which correspond to the AB 



branch of Fig. |3J and assume without loss of generality 
that the vortex is located on the x axis. We may then 
insert P = (0,P),R = (R, 0) and Q = (Q, 0) in Eqs. J2HI 
to obtain the more transparent relations 



2 Aw 

P = —, 5 Q, 

1 — UJ Z 



R = —_ 



2A 



Q. 



(29) 



The moments P, R and Q calculated along the AB branch 
(0 < I < 1) were found to satisfy the virial relations 129|) 
and are depicted as functions of angular momentum in 
Fig- El A notable fact is that all moments vanish for £ = 
or 1 because the calculated wave function becomes axi- 
ally symmetric in both of the above limits. Nevertheless, 
the moments do not vanish for other values of the angu- 
lar momentum in the interval < £ < 1. Viewed from 
the laboratory frame, the linear momentum and mean 
position of the vortex read 



-Piab = P(~ sintot, coscot), 

-Rlab = -R(C0S LOt, SlTlLUt), 



(30) 



with P — luR, in complete analogy with the motion of 
a pointlike particle rotating around the center. There- 
fore, generic behavior prevails in the presence of some 
anharmonicity (A ^= 0), whereas the stronger virial rela- 
tions (|19|l are but a curious feature of the harmonic limit 
(A = 0). 

In conclusion, the specific family of solutions analyzed 
in this paper illustrates some of the subtleties of vortex 
states in a confined Bose-Einstein condensate but cer- 
tainly does not exhaust the possibilities. It is clear that 
a huge variety of multiple vortex states are possible, with 
increasing angular momentum, which eventually lead to 
a formation of regular vortex lattices. The virial theo- 
rems derived here must be satisfied in all cases and may 
thus be used to provide important checks of consistency, 
especially because they are sensitive to the presence of 
anharmonicity. 
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